/*
 * PitzerTris.cpp
 *
 *  Created on: 2 Aug 2011
 *      Author: Allan
 */

#include "PitzerTris.h"

// GeoReact includes
#include "Geochemistry/ElectrolyteSolution.h"
#include "Utils/Algorithms.h"

const double PitzerHydron(const ElectrolyteSolution& sol, double T, double P, const VectorXd& n);
const double PitzerTRISH(const ElectrolyteSolution& sol, double T, double P, const VectorXd& n);
const double PitzerTRIS(const ElectrolyteSolution& sol, double T, double P, const VectorXd& n);

const ActivityCoefficient PitzerHydronModel(const ElectrolyteSolution& sol)
{
	return std::bind(PitzerHydron, sol, _1, _2, _3); 
}

const ActivityCoefficient PitzerTRISHModel(const ElectrolyteSolution& sol)
{ 
	return std::bind(PitzerTRISH, sol, _1, _2, _3); 
}

const ActivityCoefficient PitzerTRISModel(const ElectrolyteSolution& sol)
{ 
	return std::bind(PitzerTRIS, sol, _1, _2, _3);
}

/// The Debye-Huckel limiting law given by equation (7) 
const double fgamma(double T, double I);

/// The equivalent molality of the aqueous solution (E = 1/2 sum|zi.mi|)
const double E(const ElectrolyteSolution& sol, const VectorXd& n);

/// The parameters (B) and (C) given by equations (9) and (10)
const double B_NaCl(double T, double I);
const double C_NaCl(double T);

const double B_HCl(double T, double I);
const double C_HCl(double T);

const double B_TRISHCl(double T, double I);
const double C_TRISHCl();

/// The parameters (R) and (S) given by equations (11) and (12)
const double R(double mNa, double mCl, double T, double I);
const double S(double mNa, double mCl, double T);

/// The functions (f1) and (f2) given by equations (13) and (14)
const double f1(double I);
const double f2(double I);

/// The Pitzer parameters given by equations (15), (16) and (17)
const double beta0_NaCl(double T);
const double beta1_NaCl(double T);
const double Cphi_NaCl(double T);

/// The Pitzer parameters given by equations (18), (19) and (20)
const double beta0_HCl(double T);
const double beta1_HCl(double T);
const double Cphi_HCl(double T);

const double PitzerHydron(const ElectrolyteSolution& sol, double T, double P, const VectorXd& n)
{
	// The aqueous species
	const vector<string>& aqueousSpecies = sol.GetSpecies();
	
	// The indexes of the species H2O, Cl[-] and Na[+]
	const unsigned iH2O = IndexOf(  "H2O", aqueousSpecies);
	const unsigned iCl  = IndexOf("Cl[-]", aqueousSpecies);
	const unsigned iNa  = IndexOf("Na[+]", aqueousSpecies);
	
	// The molalities of the ions Cl[-] and Na[+]
	const double mCl = 55.508 * n[iCl]/n[iH2O];
	const double mNa = 55.508 * n[iNa]/n[iH2O];
	
	// The ionic strength of the aqueous solution
	const double I = mNa; // Assuming that mNaCl = mNa = mCl
	
	// Values taken from Table 1 from the reference
	const double Theta_HNa = 0.0289;
	const double Psi_HNaCl = 0.0002;
	
	// The natural log of the activity coefficient of H[+]
	const double ln_gH = fgamma(T, I) + 2.0*mCl * (B_HCl(T, I) + mCl*C_HCl(T)) + R(mNa, mCl, T, I) + S(mNa, mCl, T) + mNa * (2.0*Theta_HNa + mCl*Psi_HNaCl);
	
	// Return the activity coefficient of H[+]
	return std::exp(ln_gH);
}

const double PitzerTRISH(const ElectrolyteSolution& sol, double T, double P, const VectorXd& n)
{
	// The aqueous species
	const vector<string>& aqueousSpecies = sol.GetSpecies();
	
	// The indexes of the species H2O, Cl[-] and Na[+]
	const unsigned iH2O = IndexOf(  "H2O", aqueousSpecies);
	const unsigned iCl  = IndexOf("Cl[-]", aqueousSpecies);
	const unsigned iNa  = IndexOf("Na[+]", aqueousSpecies);
	
	// The molalities of the ions Cl[-] and Na[+]
	const double mCl = 55.508 * n[iCl]/n[iH2O];
	const double mNa = 55.508 * n[iNa]/n[iH2O];
	
	// The ionic strength of the aqueous solution
	const double I = mNa; // Assuming that mNaCl = mNa = mCl
	
	// The natural log of the activity coefficient of TRISH[+]
	const double ln_gTRISH = fgamma(T, I) + 2.0*mCl * (B_TRISHCl(T, I) + E(sol, n)*C_TRISHCl()) + R(mNa, mCl, T, I) + S(mNa, mCl, T);
	
	// Return the activity coefficient of TRISH[+]
	return std::exp(ln_gTRISH);
}

const double PitzerTRIS(const ElectrolyteSolution& sol, double T, double P, const VectorXd& n)
{
	// The aqueous species
	const vector<string>& aqueousSpecies = sol.GetSpecies();
	
	// The indexes of the species H2O, Cl[-] and Na[+]
	const unsigned iH2O = IndexOf(  "H2O", aqueousSpecies);
	const unsigned iCl  = IndexOf("Cl[-]", aqueousSpecies);
	
	// The molalities of the ions Cl[-] and Na[+]
	const double mCl = 55.508 * n[iCl]/n[iH2O];
	//const double mNa = 55.508 * n[iNa]/n[iH2O];
	
	// The molality of NaCl
	const double m = mCl; // Assuming that mNaCl = mNa = mCl
	
	// The natural log of the activity coefficient of TRIS
	const double ln_gTRIS = 0.04721*m + 0.000474*m*T - 0.002190*m*m - 0.0000799*m*m*T;
	
	// Return the activity coefficient of TRIS
	return std::exp(ln_gTRIS);
}

const double EqulibriumConstantTRISH(double T, double P)
{
	const double Tk = T + 273.15; // Temperature in kelvins
	
	// The natural log of the equilibrium constant of TRISH[+]
	//const double ln_KTRISH = -175.808 + 11250.187/Tk + 29.413 * std::log(Tk) - 0.0366699 * Tk;
	const double ln_KTRISH = -(-175.808 + 11250.187/Tk + 29.413 * std::log(Tk) - 0.0366699 * Tk);
	
	// Return the equilibrium constant of TRISH[+]
	return std::exp(ln_KTRISH);
}

const double fgamma(double T, double I)
{
	// The Debye-Huckel limiting slope
	const double A = 0.37674 + 5.040E-4*T + 3.30E-6*T*T;
	
	// The square root of the ionic strength
	const double sI = std::sqrt(I);
	
	return -A * (sI/(1.0 + 1.2*sI) + 2.0/1.2 * std::log(1 + 1.2*sI));
}

const double E(const ElectrolyteSolution& sol, const VectorXd& n)
{
	// The electrical charges of the ions
	const VectorXd zi = sol.GetIonicCharges();
	
	// The molalities of the ions
	const VectorXd mi = sol.IonicMolalities(n);
	
	// The equivalent molality
	double E = 0.5 * (zi.array() * mi.array()).abs().sum();
	
	return E;
}

const double B_NaCl(double T, double I)
{
	return beta0_NaCl(T) + beta1_NaCl(T) * f1(I);
}

const double C_NaCl(double T)
{
	return 0.5 * Cphi_NaCl(T);
}

const double B_HCl(double T, double I)
{
	return beta0_HCl(T) + beta1_HCl(T) * f1(I);
}

const double C_HCl(double T)
{
	return 0.5 * Cphi_HCl(T);
}

const double B_TRISHCl(double T, double I)
{
	// Values taken from Table 1 from the reference
	const double beta0_TrisHCl_25C = 0.03950;
	const double beta1_TrisHCl_25C = 0.20978;
	
	return beta0_TrisHCl_25C + beta1_TrisHCl_25C * f1(I);
}

const double C_TRISHCl()
{
	// Value taken from Table 1 from the reference
	const double Cphi_TrisHCl_25C = -0.00236;
	
	return 0.5 * Cphi_TrisHCl_25C;
}

const double R(double mNa, double mCl, double T, double I)
{
	return mNa * mCl * beta1_NaCl(T) * f2(I);
}

const double S(double mNa, double mCl, double T)
{
	return mNa * mCl * Cphi_NaCl(T)/2.0;
}

const double f1(double I)
{
	// The square root of the ionic strength
	const double sI = std::sqrt(I);
	
	return (1.0 - (1.0 + 2.0*sI) * std::exp(-2.0*sI))/(2.0*I);
}

const double f2(double I)
{
	// The square root of the ionic strength
	const double sI = std::sqrt(I);
	
	return (-1.0 + (1.0 + 2.0*sI + 2.0*I) * std::exp(-2.0*sI))/(2.0*I*I);
}

const double beta0_NaCl(double T)
{
	const double Tk = T + 273.15; // Temperature in kelvins
	
	return 14.3783204 + 0.00560767406*Tk - 422.185236/Tk - 2.51226677*std::log(Tk) - 0.00000261718135*Tk*Tk + 4.43854508/(680 - Tk) - 1.70502337/(Tk - 227);
}

const double beta1_NaCl(double T)
{
	const double Tk = T + 273.15; // Temperature in kelvins
	
	return -0.483060685 + 1.40677479E-3*Tk + 119.311989/Tk - 4.23433299/(Tk - 227);
}

const double Cphi_NaCl(double T)
{
	const double Tk = T + 273.15; // Temperature in kelvins
	
	return -0.100588714 - 1.80529413E-5*Tk + 8.61185543/Tk + 0.012488095*std::log(Tk) + 3.41172108E-9*Tk*Tk + 6.83040995E-2/(680 - Tk) + 0.293922611/(Tk - 227);
}

const double beta0_HCl(double T)
{
	const double Tk = T + 273.15; // Temperature in kelvins
	
	return 1.2859 - 0.0021197*Tk - 142.5877/Tk;
}

const double beta1_HCl(double T)
{
	const double Tk = T + 273.15; // Temperature in kelvins
	
	return -4.4474 + 8.425698E-3*Tk + 665.7882/Tk;
}

const double Cphi_HCl(double T)
{
	const double Tk = T + 273.15; // Temperature in kelvins
	
	return -0.305156 + 0.000516*Tk + 45.52154/Tk;
}
